ABSTRACT 


Predictive  modeling  is  an  important  part  of  evaluating  natural  attenuation  for 
groundwater  contaminant  plume  remediation.  This  paper  presents  a  numerical 
groundwater  contaminant  fate  and  transport  model  that  was  developed  within  a 
spreadsheet  program  for  this  purpose.  The  model  accounts  for  site  heterogeneity  and 
multiple  degradation  regions,  offering  flexibility  and  ease  of  application.  A  two- 
dimensional  finite-difference  approach  integrated  with  a  fourth-order  Runge-Kutta 
method  for  numerical  solution  was  used  as  the  mathematical  foundation.  The  model  was 
used  to  evaluate  natural  attenuation  for  removal  of  a  trichloroethylene  (TCE)  plume  from 
a  surficial  aquifer  containing  three  regions  with  distinctly  different  processes  for 
degradation  of  TCE.  Model  simulations  were  used  to  predict  how  far  the  TCE  plume 
would  migrate  within  the  aquifer.  Based  on  model  results  natural  attenuation  was  judged 
to  be  sufficient  to  prevent  migration  to  a  potential  receptor  of  the  TCE  plume. 
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Overview 


Evaluation  of  natural  attenuation  as  a  remedial  alternative  for  groundwater 
contaminant  plumes  is  accomplished  with  fate  and  transport  models.  Numerical  models 
for  this  purpose  can  easily  be  developed  within  spreadsheet  programs  that  are  available 
on  most  personal  computers.  These  programs  offer  an  easy  to  use  template  for 
development  of  numerical  models  that  are  capable  of  incorporating  site  heterogeneities 
and  multiple  contaminant  loss  mechanisms. 

The  following  paper  has  been  submitted  for  publication  in  Ground 
Water  Monitoring  and  Remediation.  The  paper  describes  the  development  and 
application  of  a  two-dimensional  numerical  contaminant  fate  and  transport  model  that 
was  formulated  in  a  spreadsheet  program.  The  model  was  specifically  developed  for 
application  to  the  Trio  Solvents  site  in  New  Brighton,  MN,  where  a  heterogeneous 
surficial  aquifer  is  contaminated  by  Trichloroethylene  (TCE). 

Alliant  Technologies,  Inc.  is  the  responsible  party  for  the  site  and  provided  funding 
for  the  natural  attenuation  study.  Conestoga-Rovers  and  Associates  is  the  consulting  firm 
tasked  with  managing  the  remediation  efforts  at  the  site.  After  five  years  of  groundwater 
extraction  efforts  (1991-1996)  failed  to  eliminate  the  contamination,  regulators  allowed 
the  extraction  pumps  to  be  turned  off  so  the  University  of  Minnesota  could  conduct  a 
natural  attenuation  study.  As  an  add-on  to  the  natural  attenuation  study,  this  modeling 
effort  was  conducted  in  order  to  offer  a  prediction  on  the  fate  of  the  TCE  plume. 
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ABSTRACT 


Predictive  modeling  is  an  important  part  of  evaluating  natural  attenuation  for 
groundwater  contaminant  plume  remediation.  This  paper  presents  a  numerical 
groundwater  contaminant  fate  and  transport  model  that  was  developed  within  a 
spreadsheet  program  for  this  purpose.  The  model  accounts  for  site  heterogeneity  and 
multiple  degradation  regions,  offering  flexibility  and  ease  of  application.  A  two- 
dimensional  finite-difference  approach  integrated  with  a  fourth-order  Runge-Kutta 
method  for  numerical  solution  was  used  as  the  mathematical  foundation.  The  model  was 
used  to  evaluate  natural  attenuation  for  removal  of  a  trichloroethylene  (TCE)  plume  from 
a  surficial  aquifer  containing  three  regions  with  distinctly  different  processes  for 
degradation  of  TCE.  Model  simulations  were  used  to  predict  how  far  the  TCE  plume 
would  migrate  within  the  aquifer.  Based  on  model  results  natural  attenuation  was  judged 
to  be  sufficient  to  prevent  migration  to  a  potential  receptor  of  the  TCE  plume. 

KEYWORDS:  natural  attenuation,  numerical  modeling,  spreadsheet,  finite-difference, 
Runge-Kutta 
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Introduction 


Natural  attenuation  is  recognized  as  a  viable  means  for  remediating  contaminated 
groundwater  (U.S.  Environmental  Protection  Agency,  1997).  The  processes  that 
constitute  natural  attenuation  include  dispersion,  adsorption,  biodegradation,  chemical 
reactions,  and/or  radioactive  decay  (Suthersan,  1997).  Although  these  processes  occur  at 
most  contaminated  sites,  they  may  be  environmentally  significant  at  only  some  (Brady  et 
al.,  1998).  Thus,  natural  attenuation  must  be  evaluated  for  effectiveness  before  it  can  be 
considered  as  a  remedial  option. 

An  appropriate  evaluation  involves  the  documentation  of  actual  contaminant  loss  at 
the  site,  as  well  as  laboratory  studies  in  which  the  mode(s)  of  attenuation  are  identified 
(Wiedemeier  T.H.  et  al.,  1996).  If  the  data  indicate  that  natural  attenuation  is  a  feasible 
remedial  alternative,  then  models  for  contaminant  transport  and  removal  are  necessary  in 
order  to  predict  the  in  situ  contaminant  fate  (Brady  et  al.,  1998).  The  predictions  can  then 
be  used  to  evaluate  the  timeliness  of  the  process  and  the  potential  for  transport  of 
contaminants  to  downgradient  receptors. 

Many  of  the  models  currently  used  for  making  these  predictions  are  limited  to 
homogeneous  sites,  including  several  that  are  based  on  the  finite-difference  method 
which  is  generally  accepted  as  the  simplest  and  most  intuitive  (Bear  and  Verruijt,  1987). 
However,  most  sites  are  not  homogeneous,  with  considerable  variation  in  the  parameters 
that  affect  transport  and  degradation.  Models  based  on  the  finite-difference  method  can 
take  these  heterogeneities  into  account  by  a  descritization  process  that  produces  an  array 
of  coupled  domains;  each  domain  is  treated  as  homogenous  with  its  own  parameters  for 
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transport  and  degradation.  Taken  together,  the  individual  domains  yield  a  migration 
profile  of  the  contaminant  plume. 

This  study  describes  the  development  and  application  of  a  two-dimensional  (2-D) 
finite-difference  spreadsheet  model  used  to  evaluate  natural  attenuation  of  groundwater 
contaminant  plumes  at  heterogeneous  sites.  The  model  was  applied  to  a  field  site  where  a 
trichloroethylene  (TCE)  plume  is  migrating  through  a  shallow  unconfined  aquifer  and 
wetland  system  towards  a  freshwater  lake. 

Model  Theory  and  Background 

In  saturated  porous  media,  contaminants  are  subject  to  advective  and  dispersive 
transport  as  well  as  sorption,  desorption  and  removal  by  degradation.  For  two- 
dimensional  transport,  the  general  equation  used  to  describe  these  influences  when 
subjected  to  unidirectional  flow  is 

dC  =  D^#C  +  D^tfC__v^dC_kC 

dt  R  dx 2  R  dz2  R  dx 

where  C  =  concentration  of  solute  in  liquid  phase  (ML'3);  t  is  time;  vx  =  average  linear 
groundwater  velocity  (LT'1);  Dx  and  Dz  are  longitudinal  and  transverse  hydrodynamic 
dispersion  coefficients  (L2T'');  R  =  retardation  factor;  k  =  first  order  reaction  constant 
(T'1)  (Javandel  et  al.,  1984).  Horizontal  and  vertical  distances  are  designated  by  x  and  z, 
respectively.  In  the  unit  designations  M  =  mass,  L  =  length,  and  T  =  time. 

The  retardation  factor  (R)  within  Equation  1  accounts  for  the  sorption  and  desorption 
of  solutes  to  aquifer  sediment;  for  aliphatic  hydrocarbons  such  as  TCE  R  can  be 
described  (Chiou  et  al.,  1979)  by  the  linear  relationship 
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r  =  1+Pj£l 

e 


(2) 


where  IQ  =  equilibrium  distribution  coefficient  (L3M'*),  pb  =  dry  bulk  mass  density  of 
the  sediment  (ML'3);  and  0  =  porosity  of  the  porous  medium  (Fetter,  1993). 

In  fate  and  transport  models  a  finite-difference  form  of  Equation  1,  obtained  by 
substituting  algebraic  differences  for  each  partial  derivative,  is  often  used  to  approximate 
the  solution  (Bear  &  Verruijt,  1987).  In  this  form  Equation  1  is 


AC,-  _  Dx  (Ci+]j  2 C/j  +  C„,y)  |  -Pz  (C,  y+1  2C,.  ;+C;  ;_!)  _  {Ci  j  CM  j)  _  ^  ^ 


At  R 


Axz 


R 


Az 


R 


Ax 


where  At  =  length  of  time-step;  Ax  =  length  of  segment  in  x  direction;  Az  =  length  of 
segment  in  z  direction.  Subscripts  i  and  j  count  cells  in  the  x-  and  z-directions, 
respectively. 

When  using  Equation  3  the  subsurface  region  (model  domain)  is  divided  into  a  grid  of 
discrete  segments  that  can  account  for  site  heterogeneity.  After  contaminant 
concentrations  along  the  domain  boundary  and  within  each  segment  at  the  beginning  of 
the  computation  (initial  condition)  are  established,  Equation  3  is  used  to  calculate  the 
change  in  concentration  within  each  segment  over  time.  The  model  domain  and 
corresponding  equations  can  easily  be  represented  within  computer  spreadsheets,  with 
each  segment  of  the  domain  depicted  by  a  spreadsheet  cell.  Separate  sheets  are 
designated  for  the  transport  and  degradation  parameters  and  each  sheet  has  a  sole 
parameter  value  for  each  segment.  The  manual  iteration  feature  of  the  spreadsheet 
program  can  be  used  to  calculate  the  change  in  concentration  within  a  cell  after  each 
time-step  based  on  concentrations  calculated  in  the  previous  time-step. 
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The  discretization  process  used  in  the  finite  difference  method  induces  numerical 
dispersion  similar  to  physical  dispersion  (Bajracharya  and  Barry,  1994).  The  Peclet 
number  and  Courant  number  are  dimensionless  ratios  that  are  used  as  criteria  for  finding 
the  optimal  segment  and  time-step  length  in  order  to  limit  numerical  dispersion  and 
maintain  model  stability  (Rao  and  Hathaway,  1 989).  The  Peclet  number  (Pe),  a  ratio  of 
the  advective  to  diffusive  terms  in  the  transport  equation,  is  defined  as 


v  Ax 
Pe  =  -^ — 


(4) 


Dx 

The  Courant  number  (Co)  is  the  ratio  of  the  advective  to  time  dependent  terms  in  the 
transport  equation  and  is  defined  as 


Co  = 


Ax 


(5) 


Generally,  numerical  dispersion  is  minimized  and  model  stability  maximized  by 
maintaining  a  Peclet  number  less  than  or  equal  to  two  and  a  Courant  number  less  than  or 
equal  to  unity  (Pinder  and  Gray,  1977).  If  these  criteria  are  met,  the  finite  difference 
model  can  nearly  result  in  the  same  solutions  for  Equation  1  as  analytical  (exact) 
solutions.  If  the  Peclet  or  Courant  numbers  are  larger  than  the  criteria,  model  results 
must  be  closely  scrutinized. 

We  eliminated  the  need  to  iterate  to  a  specific  convergence  criterion  in  each  time  step 
and  improved  the  accuracy  of  the  model’s  approximations  by  solving  Equation  3  in  each 
segment  by  a  fourth-order  Runge-Kutta  method.  The  fourth-order  Runge-Kutta  method 
was  applied  in  the  form 
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CM  =  Cl  +  \(k{  +  2  k2  +  2k3  +  k4), 

K  =hf{Ci) 

*2  =¥(<?,+!*,)  (4) 

*3=V(C,+^2) 

k4=hf(Ci+k3 ) 

where  h  =  length  of  time-step  and  f(Q)  is  the  derivative  given  by  the  right-hand  side  of 
Equation  3  (Zill,  1997). 

Model  Verification 

Output  from  the  spreadsheet  model  was  compared  to  analytical  solutions  to  establish 
model  accuracy  and  to  verify  that  the  basic  equations  were  being  solved  correctly. 
Analytical  solutions  given  by  Javandel  et  al.  (1984)  for  two-dimensional  contaminant 
transport  in  porous  media  with  one-dimensional  uniform  flow  were  used  for  the 
comparison.  The  example  consisted  of  a  unidirectional  flow  system  with  a  continuous 
conservative  contaminant  source  (Figure  1)  and  normalized  transport  parameters  that  are 
dominated  by  longitudinal  dispersion. 

Spreadsheet  model  simulations  were  run  for  four  cases  (Table  1).  All  four  cases 
closely  matched  the  analytical  solutions  given  by  Javandel  et  al.  (Figures  2  and  3).  Cases 
1  and  2  provided  the  best  fit  in  the  longitudinal  direction  because  of  the  smaller 
longitudinal  segment  length.  Cases  1  and  3  provided  the  best  fit  in  the  transverse 
direction  due  to  the  smaller  transverse  segment  length.  Case  4  seemed  to  suffer  an 
increase  in  numerical  dispersion  in  both  directions.  A  simulation  using  the  conditions 
listed  for  Case  1,  except  that  At  was  set  equal  to  0.1  day  (Co  =  0.005),  was  also 
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Table  1:  Discretization  and  Transport  Parameters  for  Verification  Problem 


Longtudnai  Distance, 


Figure  2:  Comparison  of  four  cases  of  spreadsheet  model  output  to  analytical  solutions 
adapted  from  Javandel  et  al.  (1984)  for  the  flow  domain  defined  in  Figure  1. 
Longitudinal  concentration  profiles  along  the  x-axis  from  model  simulations. 
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Transverse  Distance, 


Figure  3:  Comparison  of  four  cases  of  spreadsheet  model  output  to  analytical  solutions 
adapted  from  Javandel  et  al.  (1984)  for  the  model  domain  defined  in  Figure  1. 
Transverse  profiles  along  the  z-axis  at  x  =  50  m  from  model  simulations. 


completed  (results  not  shown).  Results  did  not  vary  noticeably  from  results  with  At  =  1. 
However,  in  simulations  completed  with  the  Courant  number  greater  than  0.2,  the  model 
became  unstable. 

Modeling  of  a  dispersion-dominated  problem  while  satisfying  the  Peclet  criterion  was 
easy,  however  if  the  example  was  advection-dominated  (i.e.  vx  =  1  md'1  and  Dx  =  0.1 
m2d_1)  the  spatial  step  size  needed  to  be  less  than  or  equal  to  0.2  m  in  order  to  meet  the 
criterion.  This  results  in  a  model  domain  consisting  of  at  least  1000  segments  along  the 
longitudinal  flow  direction  as  compared  to  10  or  more  segments  in  the  dispersion 
dominated  example.  Such  an  increased  number  of  segments  increases  model  simulation 
time. 

Model  Application  to  the  Trio  Solvents  Field  Site 

The  model  was  applied  to  a  former  solvent  distillery  site  in  New  Brighton,  Minnesota 
where  a  surficial  aquifer  is  contaminated  with  trichloroethylene.  A  groundwater 
extraction  system  was  in  operation  from  1991-1996  removing  approximately  70  kg  of 
TCE,  but  the  system  was  shut  down  after  initial  average  groundwater  TCE  concentrations 
of  2200  ppb  leveled  off  near  200  ppb.  Two  years  later,  the  average  TCE  concentration  in 
the  groundwater  was  300  ppb.  The  plume  of  residual  TCE  is  migrating  northeast  from 
the  source  area  impacting  a  stand  of  cottonwood  trees  and  then  a  wetland  that  is 
connected  to  a  freshwater  lake  (Figure  4). 

The  surficial  aquifer  is  composed  of  fine  to  medium  grained  sand  and  silty  sand  with 
a  depth  of  approximately  10  m  and  a  water  table  varying  between  0  to  1  m  below  the  land 
surface  (Delta,  1987;  Conestoga-Rovers  and  Associates,  1990).  TCE  transformation 
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Figure  4:  Aerial  View  of  Trio  Solvents  Site  showing  TCE  Source  Area,  Groundwater 
Flow  Direction,  and  Model  Domain  (A  to  A’).  (Adapted  from  Conestoga-Rovers  and 
Associates,  1990). 
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products  are  present  in  the  groundwater  (cis- 1 ,2-dichloroethylene  and  vinyl  chloride)  and 
laboratory  studies  showed  evidence  of  TCE  removal  by  indigenous  processes  (Bankston 
et  al.  1999),  indicating  that  site  conditions  are  favorable  for  natural  attenuation. 

The  actual  amount  of  TCE  spilled  in  the  past  is  unknown,  so  determination  of  current 
source  amounts  was  impossible.  Therefore  a  continuous  source  of  TCE  using 
concentration  profiles  obtained  from  groundwater  measurements  collected  over  the  last 
two  years  (Bankston  et  al.  1999)  was  used  as  input  to  the  model  domain. 

The  model  domain  was  defined  as  the  area  between  the  observed  leading  edge  of  the 
plume  and  Rush  Lake  (A  to  A’  in  Figure  4).  Bankston  (1998)  identified  three  regions 
within  this  domain  having  distinctly  different  processes  for  the  degradation  of  TCE.  A 
clay  lens  about  4  m  below  the  surface  was  also  identified.  Table  2  gives  a  summary  of 
the  estimated  ranges  for  transport  and  degradation  parameters  and  Figure  5  gives  the 
vertical  configuration  of  the  aquifer  assumed  for  the  model  domain. 

Prior  to  running  the  2-D  simulations  for  the  entire  domain,  one-dimensional 
simulations  of  a  TCE  concentration  profile  from  A  to  100  m  downgradient  were  run  for  a 
6500  day  period  and  three  different  segment  lengths  in  order  to  investigate  how  Peclet 
numbers  affected  model  output  (Figure  6).  A  simulation  that  used  a  segment  length  of 
0.012  m  to  meet  Pe  <  2  criterion  was  not  run,  because  the  number  of  segments  needed  to 
represent  the  domain  exceeded  available  desktop  computer  memory  capabilities. 
However,  other  models  have  yielded  good  results  with  Peclet  numbers  as  high  as  ten 
(Elnawawy  &  Valocchi,  1990).  Therefore  a  simulation  with  a  segment  length  of  0.028  m 
(Pe  =  5)  was  run,  requiring  over  3500  segments  to  discretize  the  domain.  A  reduction  in 
time-step  length  to  1/10  of  a  day  (Co  =  0.15)  was  necessary  in  order  to  maintain  model 
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Table  2:  Range  Estimates  for  Transport  and  Degradation 
Parameters 


Parameter 

Low 

Estimate 

Mid- 

Range 

High 

Estimate 

Reference 

Groundwater 
velocity  (md*) 

8.3  x  10 3 

4.2x10 2 

9.1  x  10 2 

Delta  (1987) 
CRA  (1990) 

Longitudinal 

dispersion 

(m2d') 

2.3  x  10"5 

1 .2  x  1 0'4 

1  x  10^ 

Bear  (1979) 

Transverse 

dispersion 

Kd') 

6.3  x  10  ® 

1.3  x  10"5 

9.4  x  10^ 

Bear  (1979) 

Aquifer  sediment 
decay  constant 
(day’) 

.00042“ 

.0077 

.04“ 

a-Appendix 

B 

b-  Howard  et 
al.  (1991) 

Cottonwood  tree 
decay  constant 
(day1) 

.00305 

.0305 

.305 

Bankston  et 
al.  (1999) 

Cattail  and 
wetland  sediment 
decay  constant 

_ l&l 

.00234 

.0234 

.234 

Bankston  et 
al.  (1999) 
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ure  5:  Vertical  configuration  of  Trio  Solvents  aquifer  from  A  to  A’  (not  drawn  to  scale)  with:  a)  the  anaerobic  aquifer  region,  b) 
cottonwood  root  region,  c)  cattail  root  region,  and  d)  the  clay  lens. 


Figure  6:  One-dimensional  simulation  results  for  segment-length  comparisons.  Results 
show  TCE  concentration  profiles  from  A  to  100  m  downgradient  for  a  simulation  time  of 
6500  days. 


stability.  As  a  result  of  the  smaller  segment  length  and  shorter  time-step  length, 
computation  time  increased.  Results  for  0.5  m  (Pe  =  83)  and  2  m  (Pe  =  334)  segment 
lengths,  which  allowed  time-step  length  to  be  set  at  1  day,  took  less  than  1/600  of  the 
computation  time  needed  for  the  Pe  =  5  simulation.  The  point  where  TCE  concentrations 
dropped  below  the  EPA’s  maximum  contaminant  level  (MCL)  of  5  ppb  was  65,  68,  and 
72  m  downgradient  for  the  simulations  with  0.028,  0.5,  and  2.0  m  segment  lengths, 
respectively.  This  was  interpreted  to  mean  that  the  TCE  plume  would  be  degraded  well 
before  reaching  Rush  Lake,  which  was  250  m  downgradient.  Literature  that  addresses 
the  evaluation  of  natural  attenuation  emphasizes  the  need  for  conservatism  in  model 
output  (Brady  et  al.,  1998;  American  Society  for  Testing  and  Materials,  1994).  Therefore 
the  time  saved  in  calculation  was  considered  to  outweigh  the  additional  numerical 
dispersion  induced  when  using  larger  segment  lengths. 

Sensitivity  analyses  of  the  transport  and  degradation  parameters  were  also  completed 
to  identify  which  parameters  had  the  greatest  influence  on  profile  results.  The  2  m 
segment,  6500  day  simulation  discussed  above  (parameter  values  in  Table  3)  was  used  as 
the  reference  profile  for  the  sensitivity  analyses.  Rate  parameters  were  varied  by  an  order 
of  magnitude  from  those  used  in  the  reference  simulation.  Variations  of  longitudinal 
dispersion  (2.5  x  10-3  and  2.5  x  10-5  m2/d)  and  transverse  dispersion  (7.5  x  10-4  and  7.5 
x  10-6  m2/d)  gave  no  appreciable  change  from  the  concentration  profile  (data  not 
shown),  indicating  that  the  system  is  advection-dominated.  However,  simulations  in 
which  the  groundwater  velocity  (4.17  x  10-1  and  4.17  x  10-3  m/d)  and  degradation  rates 
(0.0131  and  0.000131  d-1)  were  varied  by  a  factor  of  10  gave  a  large  variance  from  the 
reference  profile  (Figure  7).  As  to  be  expected,  simulations  using  a  faster  groundwater 
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Table  3:  Transport  and  Degradation  Parameters 
Used  in  Full  2-D  Simulations 


Ax 

2m 

Az 

.667  m 

At 

1  d 

vx 

0.0417  md  1 

Dx 

2.5x10J*m2d' 

Dz 

7.5  x  1 0"5  m2d‘' 

K-  aquifer 

0.00131  d1 

K-  cottonwood 

0.01  d ' 

K-  cattail 

0.0075  d ' 
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Figure  7:  Sensitivity  analyses  results  showing  influence  of  varied  groundwater  velocity 
and  degradation  rates  on  the  distance  of  plume  migration. 
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velocity  or  a  slower  rate  of  degradation  resulted  in  relatively  higher  TCE  concentrations 
being  propagated  further  downgradient  while  simulations  with  a  slower  groundwater 
velocity  or  a  faster  rate  of  degradation  resulted  in  a  reduction  of  the  plume’s  travel 
distance. 

A  simulation  for  the  vertical  (x-z)  plane  from  0  to  10  m  below  surface  and  with  the 
parameter  values  given  in  Table  3  was  completed.  This  simulation  included  all  three 
subsurface  regions  as  well  as  the  clay  lens  (Figure  5).  Groundwater  velocity  and 
dispersion  coefficients  were  assumed  to  be  equal  everywhere  except  in  the  clay  lens  that 
was  assumed  to  be  a  no-flow  zone  with  vx,  Dx,  and  Dz  equal  to  zero.  The  source 
concentration  of  TCE  was  assumed  to  be  a  uniform  350  ppb  at  x  =  0  along  the  z-axis. 
The  model  simulation  was  continued  until  the  plume  reached  the  point  at  which  the 
natural  attenuation  equaled  the  mass  flux  of  TCE  (steady-state).  This  occurred  after  6500 
or  more  days  for  this  simulation.  Figure  8  gives  the  simulation  results  for  the  predicted 
TCE  concentration  in  the  vertical  cross-section  through  the  shallow  aquifer  shown  in 
Figure  5.  The  TCE  did  not  spread  far  downgradient  in  the  shallow  depths  of  the  aquifer 
due  to  the  higher  degradation  rates  estimated  for  the  cottonwood  and  cattail  regions.  The 
MCL  of  5  ppb  was  reached  at  a  distance  of  about  17m  downgradient  of  the  source  in  the 
upper  4  to  5  m  of  the  aquifer.  In  the  deeper  part  of  the  aquifer,  i.e.  from  6  to  10  m,  the 
MCL  was  reached  approximately  72  m  downgradient  of  the  source. 

The  2-D  model  presented  and  applied  so  far  is  for  plumes  that  do  not  spread 
significantly  in  the  horizontal  direction.  This  is  typically  the  case  further  downgradient 
from  the  source  of  pollution,  especially  in  advection-dominated  situations  with  low 
transverse  dispersion.  At  the  Trio  Solvents  site  the  field  data  indicated  that  the  plume  at 
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Figure  8:  Modeling  results  of  predicted  TCE  concentrations  in  vertical  cross-section  A- 
A’  from  Figure  3  for  steady  TCE  input. 
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the  starting  point  of  the  model  domain  (point  A  in  Figure  3)  was  approximately  30  m 
wide.  Upstream  from  point  A  and  closer  to  the  suspected  TCE  source  area  field 
observations  indicate  that  the  plume  was  as  much  as  60  m  wide.  For  an  aquifer  thickness 
of  10  m  these  dimensions  justify  a  2-D  approach  for  computations  along  the  centerline  of 
the  plume. 

However,  near  the  lateral  fringes  of  the  plume  TCE  concentrations  will  naturally 
diminish.  In  order  to  obtain  an  approximate  picture  of  the  concentration  field  in  a 
horizontal  plane  through  the  bottom  part  (10  m  depth)  of  the  aquifer,  a  simulation  was 
made  with  the  model  after  exchanging  Dz  and  Az  for  Dy  and  Ay.  The  coordinate  y  is  in 
the  horizontal  transverse  direction.  The  same  numerical  value  was  used  for  Dy  and  Dz. 
In  this  simulation  the  inflow  concentrations  to  the  horizontal  model  domain  were 
specified  as  an  approximate  Gaussian  distribution  with  a  standard  deviation  of  about  7.5 
m  and  a  maximum  concentration  of  350  ppb.  The  attenuation  coefficient  for  the  aquifer 
sediment  at  the  10  m  depth  (Table  3)  was  used  for  the  simulation.  Simulation  results  are 
shown  in  Figure  9. 

Summary  and  Conclusions 

A  numerical  spreadsheet-based  model  for  a  groundwater  contaminant  plume  in  a 
heterogeneous  aquifer  was  developed  and  applied.  The  2-D  model  incorporates 
advective  and  dispersive  transport  of  a  contaminant  as  well  as  contaminant  removal  by 
sorption  and  degradation.  Heterogeneity  in  the  aquifer  and  in  these  processes  can  be 
easily  accounted  for  within  the  spreadsheets  creating  a  very  adaptable  modeling  tool  for 
evaluating  natural  attenuation  effects. 
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Figure  9:  Modeling  results  of  predicted  TCE  concentrations  in  horizontal  cross-section 
at  10  m  depth  for  steady  TCE  input. 
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The  model  output  closely  matched  analytical  solutions  indicating  that  the  governing 
equation  for  contaminant  flow  was  solved  correctly.  However,  when  applying  the  model 
to  a  field  site,  results  also  reflect  the  accuracy  of  the  input  data.  As  the  sensitivity 
analysis  demonstrated,  changes  in  flow  velocity  and  degradation  parameters  of  the 
advection-dominated  domain  alter  model  results  significantly,  whereas  changes  in 
dispersion  coefficients  did  not.  Application  of  the  model  to  the  Trio  Solvents  site  with 
input  parameters  based  on  site  knowledge  and  conservative  estimates  of  attenuation  rates 
indicated  a  high  probability  that  natural  processes  will  attenuate  the  TCE  plume  so  that  it 
will  not  impact  Rush  Lake.  Known  heterogeneity  and  multiple  degradation  regions  were 
easily  incorporated  and  model  output  showed  which  degradation  regions  had  the  greatest 
influence  on  plume  reduction.  Expansion  of  the  model  to  three-dimensions  and  inclusion 
of  additional  attenuation  processes  such  as  volatilization  could  improve  the  applicability 


for  future  use. 
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Instructions  for  Model  Use 


1.  Determine  desired  segment  length  (Ax)  and  width  (Az)  for  discretization  of  the  model 
domain  and  the  time-step  length  (At). 

2.  Open  desired  Excel  file  containing  the  model  program.  Be  sure  to  ‘click’  the  “enable 
macro”  option. 

3.  Enter  the  values  for  Ax,  Az,  and  At  in  the  appropriate  cells  on  the  sheet  labeled  Input. 
The  number  of  rows  and  columns  needed  for  the  model  application  will  automatically 
be  calculated  below. 

4.  Copy  cells  A6  through  C6  on  the  Input  sheet  down  for  the  calculated  number  of 
rows,  plus  one  more  row.  This  extra  row  insures  that  the  last  row  of  the  model 
domain  still  has  a  valid  change  in  width  value.  This  series  of  cells  allows  for 
incorporation  of  dispersion  in  the  y-direction.  Currently  it  is  set  so  that  the  cells  do 
not  change  in  width  as  we  move  down  gradient.  However,  if  a  value  for  dispersion  in 
the  y-direction  is  known  and  the  width  of  the  plume  is  known,  spreading  of  the  plume 
in  the  y-direction  can  be  included  in  calculations. 

5.  Next,  copy  cell  B6  on  the  remaining  sheets  (Co,  C,  Po,  P,  Qo,  Q,  Ro,  R,  So,  S,  V, 

Dl»  Dt»  and  k)  over  and  down  for  the  desired  number  of  columns  and  rows. 

6.  Establish  boundary  conditions.  This  is  only  necessary  on  the  Co  and  C  sheet.  An 
upgradient  source  would  be  input  by  entering  the  source  concentration  values  in  row 
5  for  the  desired  columns.  Zeros  need  not  be  entered  around  the  flow  domain  to 
signify  a  no  flux  boundary,  however  if  there  are  source  areas  along  these  boundaries, 
they  should  be  added. 

7.  Enter  initial  conditions.  This  is  done  within  the  model  domain  on  the  Co  sheet.  Enter 
the  known  concentration  values  for  each  segment  within  the  domain. 

8.  Enter  transport  and  degradation  parameters.  Enter  velocity  (sheet  V),  longitudinal 
(sheet  DL)  and  transverse  (sheet  DT)  dispersion,  and  the  first-order  degradation  rate 
(sheet  k)  parameters  for  each  segment  within  the  domain. 

9.  Press  Ctrl+Shift+A  to  set  the  counter  and  to  calculate  the  concentration  values  for  the 
first  time-step  based  on  the  initial  conditions. 


10.  Creating  the  iteration  loop:  Highlight  the  concentration  values  from  the  model 
domain  in  sheet  C,  press  Ctrl  +  C.  Go  to  sheet  Co  and  click  on  cell  B6,  from  the  Edit 
menu  select  Paste  Special  and  select  ‘Paste  Link.’  This  creates  a  loop  where  the  next 
iteration  uses  the  previous  iteration’s  concentration  values  as  initial  conditions. 
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1 1.  One  iteration  is  equal  to  one  time-step  (At).  In  order  to  run  the  model  for  a  desired 
number  of  iterations,  go  to  the  Tools  menu,  select  Options  and  select  the  Calculation 
sheet.  In  the  maximum  iterations  space  enter  the  number  of  iterations.  Click  OK. 
Then  press  F9  to  run  the  model  simulation  for  the  entered  number  of  time-steps. 
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Percent  TCE 


Figure  Bl:  Calculation  of  first-order  decay  constant  for  TCE  loss  in  cattail  experiments. 
Data  points  are  initial  TCE  present  and  TCE  remaining  in  the  microcosms  at  the  end  of 
the  experiments  conducted  by  Bankston  et  al  (1999). 
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Figure  B2:  Calculation  of  first-order  decay  constant  for  TCE  loss  in  cottonwood  tree 
experiments.  Data  points  are  initial  TCE  present  and  TCE  remaining  in  the  microcosm  at 
the  end  of  the  experiments  conducted  by  Bankston  et  al  (1999). 
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Figure  B3:  TCE  loss  in  anaerobic  aquifer  sediment  microcosms.  Symbols  are 
Experimental  microcosms  (■)  and  autoclaved  control  microcosms  (♦).  Soil  for  this 
experiment  was  collected  from  outside  of  the  TCE  plume  at  the  Trio  Solvents  site  (i.e.  in 
clean  sediment)  at  a  depth  of  8  m. 
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Figure  B4:  TCE  loss  in  anaerobic  aquifer  sediment  microcosms.  Symbols  are 
Experimental  microcosms  (■)  and  autoclaved  control  microcosms  (♦).  Soil  for  this 
experiment  was  collected  from  within  the  TCE  plume  at  the  Trio  Solvents  site  (i.e.  in 
contaminated  sediment)  at  a  depth  of  8  m. 
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Figure  B5:  Calculation  of  first-order  decay  constant  for  TCE  loss  in  “clean”  aquifer 
sediment  microcosms  by  trendline  analysis.  Only  the  data  points  collected  during  the 
first  two  weeks  of  the  experiment  were  used  for  rate  constant  determination. 
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Figure  B6:  Calculation  of  first-order  decay  constant  for  TCE  loss  in  “contaminated” 
aquifer  sediment  microcosms  by  trendline  analysis.  Only  the  data  points  collected  during 
the  first  two  weeks  of  the  experiment  were  used  for  rate  constant  determination. 


Appendix  C 
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Table  Cl:  Initial  conditions  for  vertical  cross-section  model  simulation 
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Table  C2:  Velocity  Parameters  for  vertical  cross-section  model  domain 
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Table  C4:  Transverse  dispersion  parameters  for  vertical  model  domain 
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Table  C6:  Model  output  for  100  day  simulation 
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Table  C8:  Model  output  for  6500  day  simulation 
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